﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace FastGraphClustering_13._12
{
    class Em_Algorithm
    {
        int[] Nq; // array for num of nodes in each cluster
        // double alpha_q, PI_ql, PI_qq,likelihood;

        double[] PI_qq, Nqq, alpha_q;
        int numOfClusters, numOfNodesWholeNetwork, currenteNumOfNodes, numOfNodesToAdd;
        Matrix Nql, PI_ql, Z;// matrix for num of edges between different clusters
        Network net;
        Kernel_K_Means k_means;

        public Em_Algorithm(Network currentNetwork, Kernel_K_Means currentK_Means)
        {
            this.numOfClusters = currentK_Means.numOfClusters;
            this.numOfNodesWholeNetwork = currentNetwork.numOfNodes;
            this.currenteNumOfNodes = currentK_Means.numOfNodesKernel;
            this.Nq = new int[numOfClusters];
            this.PI_qq = new double[numOfClusters];
            this.alpha_q = new double[numOfClusters];
            this.Nql = new Matrix(numOfClusters, numOfClusters);
            this.PI_ql = new Matrix(numOfClusters, numOfClusters);
            this.Z = new Matrix(numOfNodesWholeNetwork, numOfClusters);
            this.Nqq = new double[numOfClusters];
            this.net = currentNetwork;
            this.k_means = currentK_Means;
            this.numOfNodesToAdd = 0;


            for (int i = 0; i < numOfClusters; i++)
            {
                for (int j = 0; j < numOfClusters; j++)// initialize the parameters for all clusters
                {
                    Nql.Data[i, j] = 0;
                }
                Nq[i] = currentK_Means.ClusterArray[i].numOfNodes;
                Nqq[i] = 0;
                PI_qq[i] = 0;

            }

            calcNqq();


            for (int i = 0; i < numOfClusters; i++)
            {
                PI_qq[i] = (2 * Nqq[i]) / (Nq[i] * (Nq[i] - 1));// calculate PI_qq- the probability to have edge between 2 nodes at the same cluster
                alpha_q[i] = (double)Nq[i] / currenteNumOfNodes;// calculate the probability to be assignd to each cluster
            }

            calaNql();

            for (int i = 0; i < numOfClusters; i++)
                for (int j = i + 1; j < numOfClusters; j++)
                {
                    PI_ql.Data[i, j] = Nql.Data[i, j] / (Nq[i] * Nq[j]); // calculate PI_ql- the probability to have edge between 2 nodes at diffren clusters
                    PI_ql.Data[j, i] = PI_ql.Data[i, j];
                }


            for (int i = 0; i < k_means.numOfNodesKernel; i++)//initialize vector Z for each node

                for (int j = 0; j < k_means.numOfClusters; j++)
                {
                    if (j == k_means.ClusterNode[i])// k_means.ClusterNode[i]= the cluster that node i belongs to
                        Z.Data[i, j] = 1;// if node i belongs to cluster j than Z.data[i,j]=1. the probability of node i to belong to class j
                    else
                        Z.Data[i, j] = 0;//if node i dont belong to cluster j than Z.data[i,j]=0. the probability of node i to belong to class j
                }





        }

        public void calcNqq()
        {
            for (int k = 0; k < numOfClusters; k++)//itarets as number of clusters
            {
                for (int i = 0; i < k_means.ClusterArray[k].numOfNodes; i++)// iterates as number of nodes at cluster k
                    for (int j = i + 1; j < k_means.ClusterArray[k].numOfNodes; j++)
                    {
                        if (net.adjacnyMatrix.Data[k_means.ClusterArray[k].nodesInCluster[i], k_means.ClusterArray[k].nodesInCluster[j]] == 1)
                            Nqq[k]++;
                    }

            }
        }


        public void calaNql()
        {
            for (int i = 0; i < numOfClusters; i++)
                for (int j = i + 1; j < numOfClusters; j++)
                    for (int k = 0; k < k_means.ClusterArray[i].numOfNodes; k++)
                    {
                        for (int l = 0; l < k_means.ClusterArray[j].numOfNodes; l++)
                        {
                            if (net.adjacnyMatrix.Data[k_means.ClusterArray[i].nodesInCluster[k], k_means.ClusterArray[j].nodesInCluster[l]] == 1)
                                Nql[i, j]++;
                        }
                    }
        }

        public void Run_EM(int ToAdd)
        {
            int indexOfCluster;
            numOfNodesToAdd = ToAdd;
            Main_Window.calcPosFlag = 0;

            for (int i = 0; i < numOfNodesToAdd; i++)
            {
                indexOfCluster = run_Classification_Step();// index of cluster is the cluster we assign the new node to
                run_Estimation_Step(indexOfCluster);
                this.currenteNumOfNodes++;
                this.k_means.numOfNodesKernel++;// update num of nodes in kernel at Kernel K Means class
                // this.net.numOfNodesKernel_new++;//update num of nodes in kernel at Network class
            }
            if (this.net.numOfNodesKernel_old == 0)
                this.net.numOfNodesKernel_old = this.net.numOfNodesKernel_new;
            this.net.numOfNodesKernel_new += numOfNodesToAdd;
        }



        public int run_Classification_Step()
        {
            double res = 0, bestRes = 0;
            int bestCluster = 1;
            Random rand = new Random();
            double counterZi = 0;
            // choose a random cluster

            for (int i = 0; i < this.numOfClusters; i++) // count all adges between the new node with old nodes.
            {
                for (int j = 0; j < this.k_means.ClusterArray[i].numOfNodes; j++)
                {
                    if (net.adjacnyMatrix.Data[currenteNumOfNodes, this.k_means.ClusterArray[i].nodesInCluster[j]] == 1)
                        Z.Data[currenteNumOfNodes, i]++;

                }
                counterZi += Z.Data[currenteNumOfNodes, i];
            }

            for (int i = 0; i < this.numOfClusters; i++) // sum of all adges and calculate
            {
                if (counterZi != 0)
                    Z.Data[currenteNumOfNodes, i] = Z.Data[currenteNumOfNodes, i] / counterZi;
                else
                    Z.Data[currenteNumOfNodes, i] = 0;
            }
            //   Z.Data[currenteNumOfNodes, indexOfCluster] = 1;
            //    this.k_means.ClusterArray[indexOfCluster].insertNode(this.currenteNumOfNodes);  // insert the new node to cluster with maximum number of adges


            for (int i = 0; i < this.numOfClusters; i++)
            {
                //res = Math.Log10(Z.Data[currenteNumOfNodes,i])+calcSum(i);

                if (Z.Data[currenteNumOfNodes, i] != 0)
                    res = (double)Math.Log10(Z.Data[currenteNumOfNodes, i]) + calcSum(i); // alpha_q[i]= the probability to belong to cluster i
                if (Z.Data[currenteNumOfNodes, i] == 0)
                    res = calcSum(i);
                if (Math.Abs(res) > Math.Abs(bestRes))
                {
                    bestRes = res;
                    bestCluster = i;
                }
            }
            bestRes = 0;

            this.k_means.ClusterArray[bestCluster].insertNode(this.currenteNumOfNodes);// assign the new node to the best cluster
            this.k_means.ClusterNode[currenteNumOfNodes] = bestCluster;
            return bestCluster;
        }

        public double calcSum(int i)// i is the current cluster that the new node belongs to
        {
            double res = 0, pow1 = 0, pow2;
            int indexAddedNode;
            indexAddedNode = this.currenteNumOfNodes;// the index of the node we add


            for (int l = 0; l < k_means.numOfClusters; l++)
            {
                if (l != i)
                {
                    for (int j = 0; j < k_means.numOfNodesKernel; j++)
                    {
                        if (k_means.ClusterNode[j] == l) // check if Zjl=1. if node j belongs to cluster l
                        {

                            pow1 = (double)Math.Pow(PI_ql.Data[i, l], net.adjacnyMatrix.Data[indexAddedNode, j]);
                            pow2 = (double)Math.Pow((1 - PI_ql.Data[i, l]), (1 - net.adjacnyMatrix.Data[indexAddedNode, j]));
                            if (pow1 != 0 && pow2 != 0)
                                res += (double)Math.Log10(pow1 * pow2);
                        }
                    }
                }
            }


            return res;
        }
        public void run_Estimation_Step(int index)// index is the number of cluster we add the new node to
        {
            //update the Nq parameter//

            this.Nq[index]++;

            //update the Nql parameter//

            for (int i = 0; i < currenteNumOfNodes; i++)
                if ((net.adjacnyMatrix.Data[currenteNumOfNodes, i] == 1) && (k_means.ClusterNode[i] != index))
                    Nql.Data[index, k_means.ClusterNode[i]]++;

            //update the Nqq parameter//

            for (int i = 0; i < currenteNumOfNodes; i++)
                if ((net.adjacnyMatrix.Data[currenteNumOfNodes, i] == 1) && (k_means.ClusterNode[i] == index))
                    Nqq[index]++;


            //update the alpha_q parameter//

            this.alpha_q[index] = (double)(Nq[index]) / (currenteNumOfNodes + 1);

            //update PI_ql parameter//

            for (int i = 0; i < currenteNumOfNodes; i++)
                if ((net.adjacnyMatrix.Data[currenteNumOfNodes, i] == 1) && (k_means.ClusterNode[i] != index))
                    PI_ql.Data[index, k_means.ClusterNode[i]] = Nql.Data[index, k_means.ClusterNode[i]] / (Nq[index] * Nq[k_means.ClusterNode[i]]);

            //update PI_qq parameter//

            PI_qq[index] = ((2 * Nqq[index]) / ((Nq[index]) * (Nq[index] - 1)));
        }
    }
}